class BoundaryCondition {
  // supplies statevar values and possibly a volume flux or velocity profile.
  // the main model is responsible for figuring out what to do with it.
  // can be used as a true boundary condition, as a lateral river input, or
  // a mechanism for linking one Slice with another.
  // by default (until addVar() is used), does nothing. Returning null should be
  // treated by the main model as "no particular opinion: use a zero-gradient
  // condition or something equally simple."

  Slice env;
  float Qr = 0;
  float[] u = null;
  String[] names = new String[0];
  float[] values = new float[0];
  
  BoundaryCondition() {
  }
  
  BoundaryCondition(Slice env) {
    this.env = env;
  }
  
  BoundaryCondition(Slice env, float Qr) {
    this.env = env;
    this.Qr = Qr;
  }
  
  void addVar(String name, float value) {
    names = (String[])append(names,name);
    values = (float[])append(values,value);
  }
  
  int getVarIndex(String name) {
    for (int i=0; i<names.length; i++) if (names[i].equals(name)) return i;
    return -1;
  }
  
  float[] getVar(String name) {return getVar(name, env.time.current[0]);}
  float[] getVar(String name, float ti) {
    int vi = getVarIndex(name);
    if (vi==-1) return null;
    float[] result = new float[env.K];
    for (int k=0; k<env.K; k++) result[k] = values[vi];
    return result;
  }
  
  void setVar(String name, float val) {
    int vi = getVarIndex(name);
    values[vi] = val;
  }
  
  float getVolFlux() {return getVolFlux(env.time.current[0]);}
  float getVolFlux(float ti) {
    return Qr;
  }
  
  void setVolFlux(float Qr) {
    this.Qr = Qr;
  }
  
  void setVelocity(float[] u) {
    this.u = u;
  }
  
  void setVelocityFromVolFlux(int m, float h) {
    // distributes the current vol flux over the depth range -h..0.
    float Q = getVolFlux();
    float[] v = new float[env.K];
    for (int k=env.kbot[m]; k<env.K; k++) {
      float Axz = env.dx[m] * env.dz[k];
      if (-h > env.z_w[k+1]) { // too deep
        v[k] = 0;
      } else if (-h < env.z_w[k]) { // fully in input layer
        float dQ = env.dz[k] / h * Q;
        v[k] = dQ / Axz;        
      } else { // partially in input layer
        float dQ = (env.z_w[k+1] - (-h)) / h * Q;
        v[k] = dQ / Axz;
      }
    }
    setVelocity(v);
  }
  
  float[] getVelocity() {return getVelocity(env.time.current[0]);}
  float[] getVelocity(float ti) {
    return u;
  }

}



// --------------------------------------------------------------------



class NetcdfBoundaryCondition extends BoundaryCondition {
  
  // reads a netcdf file and saves time-depth fields for some number of statevars
  
  float[][][] tzvalues;
  float[] tbc, zbc;
  float[] volFlux = null;
  
  NetcdfBoundaryCondition() {}
  
  NetcdfBoundaryCondition(Slice env, String ncfilename) {
    this.env = env;
    NetcdfFile nc = nc_open(ncfilename);
    tbc = nc_read1D(nc,"t");
    zbc = nc_read1D(nc,"z");
    tzvalues = new float[0][tbc.length][zbc.length];
    // look for statevars in the netcdf file
    for (int vi=0; vi<env.vars.length; vi++) {
      String varname = env.vars[vi].shortname;
      Variable v = nc_getVar(nc, varname);
      if (v != null && !varname.equals("t")) { // found a statevar
        float[][] data = nc_read2D(nc, varname);
        tzvalues = (float[][][])append(tzvalues, data);
        names = (String[])append(names, varname);
      }
    }
    // look for volFlux in the netcdf file
    Variable v = nc_getVar(nc, "volFlux");
    if (v != null) volFlux = nc_read1D(nc, "volFlux");
    nc_close(nc);
  }
  
  float[] getVar(String name, float ti) {
    int vi = getVarIndex(name);
    if (vi==-1) return null; // not found so return null, as if the b.c. weren't defined
    // having found it...
    float[][] data = tzvalues[vi];
    float[] result = new float[env.K];
    ti = constrain(ti, tbc[0], tbc[tbc.length-1]); // constraining ti and zi is equivalent to extrapolating the first and last b.c. values  to points beyond the range where the b.c. is defined
    for (int k=0; k<env.K; k++) {
      // at each grid level, interpolate into the bcdata
      float zi = constrain(env.z_rho[k], min(zbc[0],zbc[zbc.length-1]), max(zbc[0],zbc[zbc.length-1]));
      result[k] = interp2(tbc, zbc, data, ti, zi);
    }
    return result;
  }
  
  float getVolFlux(float ti) {
    if (volFlux==null) return 0;
    ti = constrain(ti, tbc[0], tbc[tbc.length-1]);
    return interp1(tbc, volFlux, ti);
  }
  
  // if u arrived as a statevar from the netcdf file, use that. If not, use the special internal store 
  float[] getVelocity(float ti) {
    float[] uu = getVar("u",ti);
    if (uu==null) {
      return super.getVelocity(ti);
    } else {
      return uu;
    }
  }
  
}

